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We generalize the uniform-gas correlation energy formal- 
ism of Singwi, Tosi, Land and Sjolander to the case of an 
arbitrary inhomogeneous many-particle system. For jellium 
slabs of finite thickness with a self-consistent LDA ground- 
state Kohn-Sham potential as input, our numerical results for 
the correlation energy agree well with diffusion Monte Carlo 
results. For a helium atom we also obtain a good correlation 
energy. 
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Despite eighty years of study, the accurate calculation 
of the correlation energy of interacting quantal many- 
electron systems is still often a challenge, even for some 
systems not regarded as "strongly correlated" . For real- 
istic many-electrons systems the current state of the art 
includes diffusion/ Green function quantum Monte Carlo 
(DMC), variational and quantum chemical methods such 
as the Configuration Interaction (CI) approach. These 
all have stringent practical limitations to relatively small 
and/or not-too-complex systems. Approximate density 
functionals of the local-density (LDA) and generalized 
gradient (GGA) class 0] are in principle less accurate 
than the above approaches, but they remain feasible even 
for very large, complex systems, and often provide use- 
ful accuracy. They fail completely, however, to describe 
long-ranged correlations in cases where these differ qual- 
itatively from those of the homogeneous electron gas. A 
case in point is the van der Waals (vdW) or dispersion in- 
teraction: all LDA/GGA approaches miss its long-ranged 
part completely 0] JI Q, and give at best patchy results 
at short range [§, [gfT Thus DMC, CI and standard DFT 
methods are all likely to be problematical for large com- 
plex vdW systems of practical interest, including many 
soft-matter cases. 

Here we present a new general approximation method 
("ISTLS") for the correlation energy of inhomogeneous 
electronic systems, which we believe will be appropri- 
ate for vdW problems among others. It employs a self- 
consistent scheme for the pair distribution based on the 
surprisingly accurate homogeneous-gas correlation energy 
method of Singwi, Tosi, Land and Sjolander (STLS) 0. 
We take as input an approximate Kohn-Sham (KS) po- 
tential vks {f) of the inhomogeneous system, and produce 
the correlation energy E c as output. Tractable input the- 
ories such as LDA or exchange-only Krieger-Li-Iafrate 
, by themselves, make large correlation-energy errors 



of order +100% and -100% respectively, yet starting from 
their KS potentials our ISTLS theory yields correlation 
energies accurate to a few percent. In its numerical com- 
plexity, and also in its accuracy for E c of large systems, 
ISTLS appears to be intermediate between LDA/GGA 
and the more microscopic approaches mentioned above. 

Our original motivation for deriving ISTLS was pri- 
marily to address soft-matter problems such as polymer 
cohesion or the energetics of graphite and its intercalates. 
For such vdW systems, the RPA-like || nature of the 
method, together with its uniquely self-consistent local- 
field correction, suggests it will accurately describe vdW 
interactions J2| while also correctly treating other kinds 
of bonding [|4|. The method is by no means restricted 
to vdW problems, however, and should provide a useful 
alternative approach both for soft and for hard matter. 
In particular it appears to be highly competitive with 
recently-proposed correlation energy theories based on 
the GW method || ; it gives a much more accurate cor- 
relation energy in the homgeneous gas at large r s , for ex- 
ample. Being intrinsically approximate, ISTLS needs to 
be tested. Therefore, in the present paper, we benchmark 
our method against state-of-the-art results for two sim- 
ple but highly inhomogeneous situations: finite-thickness 
jellium slabs and the He atom. The results suggest that 
our scheme, while based on approximations known to 
work well in uniform electron gases, also gives a good 
treatment of strong inhomogeneity in one to three space 
dimensions. Details follow. 

The groundstate energy of an inhomogeneous many 
electron system with external potential v ext (r) and 
groundstate electron density n(f) is given exactly by the 
constant-density adiabatic connection formula of formal 
Kohn-Sham density functional theory: [p^| , p| 

E = T s [n] + J n(r)v ext (r)dr (1) 

1 r e 2 

+ - / F ^M r)n(f')drd? + E xc [n] (2) 

2 J \r - r'\ 

If 1 f e 2 
E xc [n] = - d\ — — — (n 2 \(f, r) - n(r)n(r'))drdr'. 
2 Jo J \r-r'\ 

(3) 

Here T s [n] = h 2 {2m)' 1 Y, kocc f |V0 fc (r)| 2 dr is the KS 
kinetic energy, and {</>fc}are the occupied KS orbitals, 
eigenfunctions of the one-electron KS potential Wifs(r). 
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vks is defined to be such that independent electrons mov- 
ing in vks yield the true groundstate density: 



v K s => {4>k(r)} ■ n(r) = 



k occ 



(4) 



The A integration in (J3J) accounts for the kinetic part of 
the KS correlation energy. The groundstate pair distribu- 
tion n,2\(r, r ') is that of a "A- system "defined to have a re- 
duced electron-electron interaction Ae 2 /ri2, and a modi- 
fied external potential v x xt (r) chosen to maintain the true 
(A = 1) groundstate density at any A: n x (r) = n x =\{f) = 
n(r). Remarkably |l0| ], |11|, only the true external po- 
tential v ext = v^ x appears in (0). Note that, from (Q), 
the KS potential of each A-system is the same as that of 
the true (A = 1) system, because the density is the same: 



vks,\ = vks- 



(5) 



The groundstate pair distribution n 2 \ in (||) can be re- 
lated to the Kubo density-density response function xx of 
the A-system by the T = OK fluctuation-dissipation the- 
orem @, @, 

n(r)n(r')(g\(r,r") - 1) = n 2X (r 1 r") - n(f)n(r") 

c oo 

X\(r*,r ,oj = iu)du — n(r)5(f — r* ). (6) 

Eq. (^|) also introduces the equilibrium pair correlation 
factor g x (r, r '). 

The Kohn-Sham density-density response \ks @ ^ s 
defined to be that of independent electrons moving in the 
KS potential vks- Note that, by (||), 



Xks,x = Xks — Xx=o- 



(7) 



Xks is exactly expressible |14| by perturbation theory in 
terms of the KS orbitals {4>k}- From this expression it 
is readily shown that, when xx is replaced by xa=o = 
Xks, Eq. (||) gives the "exact DFT exchange" energy, 
i.e. it gives the Hartree-Fock energy integral in which 
the self-consistent Hartree-Fock orbitals are replaced by 
the KS orbitals cf>k ■ Thus this formalism easily deals with 
exchange. Subtracting this DFT exchange energy from 
(0) we obtain the exact DFT correlation energy 



-h f . e 2 
E c = tt" / drdf'— — - 
2n I r — r 



|15[ IQ. Here, however, instead of using a local uniform- 
gas-based approximation we effectively generate a nolo- 
cal f xc self-consistently for the particular inhomogeneous 
system. To do this we extend to nonuniform systems 
the semiclassical approach of STLS. Thus we relate the 
independent-electron and interacting responses by solv- 
ing the time evolution equation (first BBGKY hierarchy 
equation jlTl) for the one-electron distribution function 
f(f,p,t) of the classical A-system: 

<I=F^> ' i-J?'^.V-,?',v',t)dr'diS'. (9) 
or \r — r'\ op 

This equation is exact but requires the dynamic pair dis- 
tribution f{ '. The essential contribution of STLS was 
to use the equilibrium pair-density factor g x (r, r') of Eq. 
(^) in a semi-classical truncation scheme 

f x 2) {r,p;r',p',t) « gx(r,f')f(r,p,t)f(f',p',t), (10) 

where the true dynamic correlation factor g\ should de- 
pend on both the momenta and the time, but this de- 
pendence is ignored and g\ is taken to be the static, 
momentum-independent equilibrium density correlating 
factor from Eq. (^|). Using ( |l0| ) in (^) and lineariz- 
ing about the equilibrium distribution, / = fo(r,p) + 
Sf(r,p, t), we obtain a closed one-body kinetic equation 



Sf(r,P,t) 



-6F eff (r,t) 



dfo(r,p) 
dp 



(11) 



Here 



F (0) (r) = F T(r) - J (^7^77 )9xfr ? >o(r ')dr' 

corresponds to the gradient of the KS potential in the 
quantal case, and is independent of A by choice of (r). 
Further, 



8F eff (r,t) = SF ext (r,t) + / W x (r,f')Sn(r' ,t)dr' , (12) 

(13) 



W x {f,r')=g x {T,r')- ° 



x d\ du[x\(r 7 r',uj = iu) - xx=o{r,r' ,u = iu)}. 
Jo Jo 

At each A value the interacting and KS response are re- 
lated exactly by a Dyson-like screening integral equation 
@ Xx = XA=o + XA=o*(AV co „/ om6 + /z C A)*XA,where spa- 
tial convolution is represented by a star. The xc kernel 
fxcX contains the many-body xc effects and has tradion- 
ally been treated by a local density approximation O] 



dr \r — r'\ 

(8) Because (O) is linear and time invariant, its solution 
Sn(r,t) = j6f(r,p,i)dp can be expressed in the form 

Sn(r,t) = I ' v {r,r l 1 t~t l ).5F en {r' 1 t')dr"dt l (14) 



where a z?o is a A-independent classical vector response 
function giving the independent-electron density re- 
sponse to an applied force, with (d/dr').i>o — xa=o- 
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In the case of homogeneous electron gases, Fq is zero 
and g\ is a function only of the separation R = \r — r'\, 
and then V x W = so that the effective pair force 
W\ is irrotational and can be expressed as a gradient 
of a scalar potential, W\{R) — —(d/dR)w\{R). Then, 
assuming 5F ext comes from a potential SV ext we can 
use integration by parts (Green's theorem) followed by 
space Fourier transformation to write ( |l4| ) in q— space as 
8n = x\=o5V eff = x\=o(SV ext + w x Sn). This yields 
X\(q,v) = XA=o(?,w)(l - w x {q)x\=o{q,^))^ 1 - This 
equation resembles a classical Random Phase Approxi- 
mation (RPA), with W\(q) replacing the bare coulomb 
pair potential Air\e 2 /q 2 . This response xa depends on 
gx via w\ , and g\ is determined by \\ y i a @ > giving a 
closed selfconsistent scheme. STLS applied this theory 
to the degenerate electron gas by replacing the classi- 
cal Boltzmann-cquation density response xa=o with the 
quantal Lindhard response. Despite the crudeness of 
the factorization (|l0|), the STLS formalism gives excel- 
lent correlation energies for both 3D and 2D homoge- 
neous electron gases, up to relatively large values of the 
inter-electron spacing parameter r s . For example, in 3D 
E^ TLS is within about 1% of the 3D diffusion Monte 
Carlo results [H for 2 < r s < 5. The error is still under 
4% at r s = 20 and under 7% at r s = 50, a regime includ- 
ing gases generally regarded as significantly correlated. 
These results are significantly better than, e.g. recent 
GW-based many-body methods || which give 20% error 
at r s = 20. The homogeneous STLS scheme has a num- 
ber of shortcomings including unphysical negative values 
of the on-top pair factor g(r, r) and failure to satisfy the 
compressibility sum rule. Further work addressing these 
difficulties @, go|, g|, @, || did not, however, sig- 
nificantly improve the predicted uniform egas correlation 
energy. Therefore in the present work we have concen- 
trated on generalizing the original semi-classical STLS 
scheme to inhomogeneous systems. This does not ap- 
pear to have been attempted previously: bilayered elec- 
tron gases have certainly been treated Q, but these 
are isomorphic to a two-species homogeneous 2D elec- 
tron gas. We will show that the formalism is tractable 
for cases of genuine inhomogeneity. 

In an inhomogeneous system we haveV x W\ = 
—Vg\(r,r') x V(Ae 2 / \r - r'j) / so that there is no 
scalar potential correspon ding to W\, and the vector 
bare response vq from Eq. ( jl4| ) must be used: the scalar 
version \o is n °t sufficient. This is an essential differ- 
ence between the inhomogeneous case and the homoge- 
neous one. As in the homogeneous case, we postulate 
that a degenerate Fermi system can be treated via the 
above semi-classical analysis by using the quantal-Fermi 
independent-electron response for : this also amounts 
to using the quantal KS potential in place of its clas- 
sical counterpart Fq. By perturbation of the occupied 
independent-electron (Kohn-Sham, KS) orbitals (f>j(r) , 



we obtained for the inhomogeneous quantal response at 
imaginary frequency iu 



u (r, f', u = iu) = - Re[— V fjtf (r) 
u m 



j 



x[G(r,r',E= huj + ihu)V '<^(r') 
-<t>j(r')V'G(r,r',E= huj^+ihu)] (15) 

where fj is the Fermi occupation factor and G is the 
Green function for a single electron moving in the 
groundstate Kohn-Sham potential vks{t)- The Coulomb 
screening conditions (|l2|), (|l3|), ( |l4| ) for the inhomoge- 
neous case can be written as a Dyson-like "screening" 
integral equation for the interacting response x\'- 



X\(r,r l ,uj) = XKs{r,r',uj) + J Q x (r, r" , u))xx{f" , f , iv)dr" , 

(16) 

Q A (fy» = J ^(f/V^^'.V. Xks = V'-vo 

(17) 

We term this the "inhomogeneous STLS" (ISTLS) 
scheme. To demonstrate its feasibility and accuracy we 
have carried it out numerically for two highly inhomoge- 
neous but spatially symmetric cases, namely (i) charge- 
neutral jcllium slabs and (ii) a helium atom. 
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FIG. 1. Correlation energy of neutral jellium slabs 
(Hartree/el) 

The jellium slabs were first solved in the LDA 
(LDA-Perdew-Wang 1991 ^5|) to give the selfconsistent 
groundstate Kohn-Sham potential vj^g A (z) and density 
n(z), where z is the space coordinate in the thin dimen- 
sion of the slab. Our ISTLS formalism was applied as 
a "post-functional" giving the correlation energy start- 
ing from the fixed v^g A (z) , though of course ideally one 
would choose ai>(z)to minimize the total energy includ- 
ing the STLS corrections. (This OPM method will also 
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give an improved density n(r)). Eqs. (0), (|T^) and jl5| ) 

- ( p"7| ) were then implemented. To perform the iterative 
refinement of the static pair factor g\(r,r') , we started 
with a g of Hartree-Fock form based on the occupied KS 
orbitals from the LDA groundstate. Finally (||) was used. 

Fig. 1 gives the slab correlation energy per electron 
(open symbols) for a number of positive background den- 
sities no+, parametrized by the dimensionlcss intcrclcc- 
tron spacing r s — me 2 ?i _2 (3no+/47r) 1 / 3 . We show results 
from ISTLS (solid line), DMC §6| (dotted line) and RPA 
(dashed line) schemes, the last obtained by setting the 
pair correlation factor g to unity in Eqs. ( p"3| ) and Jl5| ) 

- (0) . The thickness of the positive background in each 
slab is L = 7.21r s aB, to match the available DMC re- 
sults. Results per electron in the unbounded uniform gas 
are also shown (closed symbols). Agreement of E^ STLS 
with the slab DMC data is good, within 3% . This is 
comparable to the agreement of STLS with DMC for the 
uniform 3D gas with 2 < r s < 20 . 

For small finite systems such as atoms one needs a 
self-interaction correction (SIC) in the starting KS po- 
tential and density. Otherwise (as for example when 
one uses the simple LDA) unrealistic response functions 
are obtained because the asymptotic —e 2 /r potential is 
missing in v KS . We solved helium using the Krieger- 
Li-Iafrate exchange-only description ||] of the atomic 
groundstate. This has the advantage of a common poten- 
tial v KLI = v KS for all orbitals. The explicit spherical 
form of Eqs. ( |l5| ) - ( |j~7| ) involves spherical harmonics but 
is somewhat cumbersome because of the vector character 
of z7q. We obtained a total ISTLS He correlation energy 
of —40.0 milliHartree, within 5% of the "exact" nonrela- 
tivistic value [^8| of —42.0 mH. Our result is of "chemical 
accuracy" . The KLI starting potential is adequate: we 
re-ran our method starting from the numerically exact 
He KS potential (^tJ, obtaining << 1% change in E c . 

In summary, we have derived an inhomogeneous gen- 
eralization ("ISTLS") of the rather successful STLS 
uniform-gas correlation energy formalism: see Eqs. (13), 
T5|), (iJ), (|l|), (|), and (|). We have shown that 
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ISTLS gives good groundstate correlation energies in 
some highly inhomogeneous electronic systems. The 
scheme can also encompass finite temperatures, and plas- 
mon calculations. An advantage of the ISTLS scheme is 
that its pair correlation physics is "self-tailored" to the 
system at hand, rather than "stolen" from some reference 
system as in some density functional schemes, or guessed 
then optimized as in variational schemes. It does this, 
moreover, without the need for solving explicit two-body 
dynamic equations. We speculate that interesting den- 
sity functionals for the correlation energy might be de- 
rived by using the ISTLS scheme with semi-local-density 
approximations for the bare response vq , somewhat as in 
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